c
c    To be linked with wp_xi_cf_hod3.f
c    This module gives the 1-halo and two halo term power spectrum
c    N_cen(M):   Erf form: 0.5*(1+erf((log M-log Mc)/Sig Log M))
c    N_sat(M):   PL form:  ((M-M0)/M_s)^alpha 
c      
      subroutine init_nm(argc2)      
      ! Read model parameters/ In this model
      ! Mcrit_1,alpha1,filename
      parameter (nmass=50,idv=70)
      character*(*) argc2
      character nm_fil*80, rline*132, cut_space*132
      real lgmcc,sig_lgmc,logms0,lgms1,alpha
      common/wgtpar/lgmcc,sig_lgmc,logms0,lgms1,alpha
c
c     Read Table for population 2
c
      read (argc2,*) nm_fil
      open (unit=70,file=nm_fil,status='old')
      do while (.true.)
         read (70,'(A)',end=3) rline
         rline = cut_space(rline)
         nm = 0
         if (index(rline,'!').ne.1) then
            nm = nm+1
            read (rline,*) lgrms(nm),gn_sat(nm),gn_cen(nm)
            if (lgrms(nm).gt.1e5) lgrms(nm)=log10(lgrms(nm))   
!             If entry is lin, covert to log
         end if
      end do        
 3    close (70)
      ! Output Header
      write(6,'(A)')
     $     '! ACF model acf_tbl/From Table'
      write(6,'(A,A,1PE9.2,0PF5.2)')
     $     '!   N(M) file:',nm_fil(1:index(nm_fil,' ')-1)
      end
c
c
      real function galn_halo (rm,isam,icen) 
      ! For the ACF isam=1 and isam =2 gives the same
      ! outputs.
      ! rm: Mass of DMH in solar mass
      ! isam: sample number
      ! icen=0: satellite galaxy, icen=1: central galaxy 
      parameter (nmass=50)
      real lgmcc,sig_lgmc,logms0,lgms1,alpha
      common/wgtpar/lgmcc,sig_lgmc,logms0,lgms1,alpha
      lgrm = log10(rm)
      if (lgrm.lt.lgrms(1)) then
         galn_halo = 0
      else
         !
         if (icen.eq.0) then    ! Sate
            galn_halo =  (1.-ff)*gn_sat(imin)+
     $           ff*gn_sat(imax)
         else                   ! Center
            galn_halo = (1.-ff)*gn_cen(imin)+
     $           ff*gn_cen(imax)
         end if
      end if
      return
      end
c
c
      real function wgt_cf(rm,icen) ! Calculate Weight for
                              ! Expected number of pairs 
                              ! between sample 1 and sample 2
      rnmax_subp = 3.0
      prob_max = 0.999 
      rm10 = rm/1e10          ! rm in units of 1e10 solar mass/to match with wgt_samp   
      if (icen.eq.0) then     ! Sat vs Sat, subpoisson behavior 
         rngal = galn_halo(rm,1,0)
         if (rngal.le.rnmax_subp) then
            do while (p_sum.lt.prob_max)
               prob = pois (rngal,nn) ! Poisson Distr
               p_sum=p_sum+prob
               rnpair = rnpair+prob*nn*(nn-1)
               nn = nn+1
            end do
         else
            rnpair = rngal*(rngal-1)
         end if
         wgt_cf =  rnpair/rm10/rm10
      else   ! Central vs satellite
         wgt_cf = (galn_halo(rm,1,1)*galn_halo(rm,1,0))/rm10/rm10
      end if   
      return
      end
c
      subroutine model_help()
      print '(A)',' ACF Model Specification: acf_tbl.f'
      print '(A)',' N(M)_pars: 5-parameter model after Zheng+05 ApJ 633,791'
      print '(A)','      LgMcc, SigLgM,LgMs0,LgMs1,alpha'
      print '(A)','   N(M)_cen = 0.5*erf((Log M-LgMcc)/SigLgM))'
      print '(A)','   N(M)_sat = ((M-Ms0)/Ms1)^alpha'
      print '(A)', ''
      return
      end
